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Abstract — There are limitations on the extent to which man¬ 
ually constructed mathematical models can capture relevant 
aspects of legged locomotion. Even simple models for basic 
behaviours such as running involve non-integrable dynamics, 
requiring the use of possibly inaccurate approximations in the 
design of model-based controllers. In this study, we show how 
data-driven frequency domain system identification methods can 
be used to obtain input-output characteristics for a class of dy¬ 
namical systems around their limit cycles, with hybrid structural 
properties similar to those observed in legged locomotion systems. 
Under certain assumptions, we can approximate hybrid dynamics 
of such systems around their limit cycle as a piecewise smooth 
linear time periodic system (LTP), further approximated as a 
time-periodic, piecewise LTI system to reduce parametric degrees 
of freedom in the identification process. In this paper, we use a 
simple one-dimensional hybrid model in which a limit-cycle is 
induced through the actions of a linear actuator to illustrate 
the details of our method. We first derive theoretical harmonic 
transfer functions of our example model. We then excite the 
model with small chirp signals to introduce perturbations around 
its limit-cycle and present systematic identification results to 
estimate the harmonic transfer functions for this model. Com¬ 
parison between the data-driven HTF model and its theoretical 
prediction illustrates the potential effectiveness of such empirical 
identification methods in legged locomotion. 

I. Introduction 

Legged locomotion emerges from a staggering diversity 
of animal and robot morphologies and gaits, and modeling 
locomotor dynamics remains a grand challenge in both biology 
and robotics [1,2]. Running behaviors, in particular, are com¬ 
monly represented by relatively simple spring-mass models 
such as the Spring-Loaded Inverted Pendulum (SLIP) model 
[3,4]. A common feature of such models, however, is that 
their hybrid system dynamics involve intermittent foot contact 
with the ground, alternating between flight and stance phases 
during locomotion. Despite the presence of seemingly simple 
models for basic behaviors such as running and walking, the 
hybrid dynamics associated with these behaviors can be rather 
complex, with non-integrable parts such as the stance phase [3, 
5], Given the utility of having accurate models and associated 
analytic solutions in constructing high performance controllers 
for nonlinear systems, substantial effort has been devoted to the 


construction of approximate solutions to such non-integrable 
hybrid models [6-9]. 

When accurate analytical solutions to the dynamics of 
a legged platform are available [8], their structure can be 
exploited to yield effective solutions for system identifica¬ 
tion and adaptive control [10]. Despite our previous studies 
showing how accurate such models may be, there will always 
be unmodeled components in the physical system, resulting 
in discrepancies between the model and experiments [11], 
Attempts to manually incorporate these effects into the model 
is daunting at best, and often impossible. Consequently, we 
propose an alternative method in this study, namely using data- 
driven system identification methods to derive an input-output 
transfer function for such hybrid legged locomotion behaviors, 
thereby eliminating the need to manually construct an explicit 
mathematical model for the system. 

Our main goal in this study is to provide a system 
identification framework applicable to a useful (although not 
comprehensive) class of legged locomotion models [8], and 
possibly more complex robotic systems [12]. Our approach 
is based on considering legged locomotion as a hybrid non¬ 
linear dynamical system with a stable periodic orbit (limit- 
cycle), corresponding to the locomotor behavior of interest. 
We introduce a formulation that addresses the input-output 
system identification problem in the frequency domain for a 
sub-class of hybrid legged locomotion models. More specifi¬ 
cally, following certain assumptions on the hybrid dynamics 
of legged systems, we approximate their hybrid dynamics 
around the limit-cycle as a linear time-periodic system (LTP). 
However, this first LTP approximation is infinite dimensional, 
making parametric identification challenging. We hence further 
approximate the dynamics as a finite dimensional piecewise 
LTI system (maintaining its LTP nature), thereby limiting 
the parametric degrees of freedom while enabling a practical 
identification framework. 

Existing studies on system identification of LTP systems 
focus on modeling these systems as multi-input single-output 
LTI systems. This approach is based on the concept of har¬ 
monic transfer functions [13], which are infinite-dimensional 
operators that are analogous to frequency response functions 



for LTI systems. An identification strategy for such systems 
was developed in [14] using power spectral density and cross 
spectral density functions. A similar method was used in [15] 
considering the effects of noise in both input and output 
measurements. Different than these studies, local polynomial 
methods and lifting approaches were also used for the identi¬ 
fication of harmonic transfer functions for multi-input single¬ 
output models of LTP systems [16]. 

Our contributions in this paper focus on representing the 
dynamics of legged locomotion as a linear time periodic 
system, thereby enabling the use of the system identification 
method proposed in [14] for such systems. We achieve this 
by using a new phase definition in identifying the harmonic 
transfer functions, illustrated in the context of a simplified 
model designed to mirror structural properties of legged lo¬ 
comotion models. When the problem is approached as a grey- 
box model with finite parameters (piecewise LTI), it suffices 
to non-parame trie ally estimate a finite number of harmonics, 
to which we later fit parametric models. 

Ankarali and Cowan [17] developed a similar system 
identification method for hybrid systems with periodic orbits 
using “discrete time” harmonic transfer functions. However, 
the framework and assumptions in this paper are distinctly 
different from their approach. Specifically, they use mappings 
between different cross sections to construct a discrete-time 
LTP system, and use discrete time HTFs for identification. 
Also, the current paper focuses on harmonic balance, which 
also distinguishes this paper from [17]. 

II. Representation of Legged Locomotion as a 
Hybrid Dynamical System 

Our goal in this study is to provide a system identification 
framework for a class of models related to legged locomotion 
using harmonic transfer functions. For the present paper, 
we limit ourselves to “clock-driven” locomotion models as 
described in Section II-A, representative of controllers used 
by a wide variety of existing robots [12,18,19], with open- 
loop central pattern generators (CPG) coordinating control 
actions to achieve time periodic behaviour. This will allow 
us to directly use time periodicity in our LTP analysis, while 
eliminating a variety of complications associated with estimat¬ 
ing the phase [20]. 


q{t ) = q(t — T) when u(t) = 0. For such systems, if we let 
q(t) = q(t) + x(t ) and linearize the dynamics in (1) around 
the limit-cycle q{t), and u{t) = 0 we get 

x(t) = A(t)x(t) + B{t)u{t) 
y{t) = C(t)x(t) + D(t)u[t) 


where 
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u{t ) — 0 
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du q(t) = q(t) 

u(t ) = 0 


(3) 

(4) 


This corresponds to a Linear Time Periodic (LTP) system, with 
all system matrices sharing a common period, T. 


B. Modeling Framework for Hybrid Systems 

Legged systems are often modeled using hybrid dynamics 
due to intermittent foot contact with the ground, which cannot 
be represented with a single, smooth dynamical flow. In 
the broadest sense, a hybrid dynamical system is a set of 
smooth flows together with discrete transitions (and associated 
transformations) between these flows triggered by intersections 
of system trajectories with sub-manifolds of the continuous 
state space [21]. These flows are called charts , indexed with 
unique labels I := {0, • • • , d} each with possibly different 
equations of motion. Along its trajectories, a hybrid system 
transitions from one chart to another, with each transition 
defined by the zero crossing of a threshold function. For each 
source chart «el and destination chart (3 £ 1, the threshold 
function h £ defines the transition from chart a to chart /3. 
An example transition graph for a hybrid dynamical system is 
illustrated in Fig. 1. 

Since we are interested in the local behaviour around 
the limit-cycle, we assume that there is only one transition 
function associated with each chart. 1 We further assume that 
system trajectories are continuous at transitions, meaning that 
system states do not experience discrete changes coincident 
with chart transitions. As a final note, we assume that the 
hybrid dynamical system we consider has an isolated periodic 
orbit ensuring that chart transitions within the limit cycle are 
also periodic and consistent. 


A. Smooth Clock-driven Oscillators 

In general, the dynamics of smooth, clock-driven oscillators 
with external inputs can be written as 

q = f(qf,u), 

<j>=l 

(1) 

/ : I" x S 1 x l p 4 R n 
(q, <f>) £ R" x S 1 , u£WL q 

where (q, (f> ) and R" x S 1 denote the state vector and the 
state space of the oscillator, respectively. The circle component 
S 1 = mod(K + ,T) enforces the periodicity of the dynamics, 
while the external input u{t) represents small external pertur¬ 
bations which we will use for system identification. 

In this paper, we focus on oscillators of the form (1) 
with asymptotically stable, isolated periodic orbits (limit-cycle) 


It is important to note that these assumptions are generally 
satisfied by models of common locomotory behaviors such as 
running and walking [8,22] as well as a wide range of legged 
robots for which leg masses are negligible compared to the 

'This approach does not apply to gaits such as pronking that nominally 
involve multiple legs making contact with the ground at the same time when 
on the limit cycle, because small deviations from the limit cycle can lead to 
different touch-down order between legs, violating our assumption. 


Fig. 1. A simple state transition graph for a hybrid dynamical system. 



dynamics of a larger body [12,18]. Consequently, the system 
identification methods we introduce will remain applicable to 
systems other than the simplified example we will present in 
this paper. 


C. Modeling Legged Locomotion as a Linear Time Periodic 
System 

For clarity, we limit our focus in this section to an example 
hybrid dynamical system with only two charts, I = {0,1}, 
designed to capture stance and flight phases of simple spring- 
mass models of locomotion. Based on a clock driven as¬ 
sumption, for each i £ I the continuous dynamics can be 
represented with 

0=1 

ii = , (5) 

ft G R" 


and let the associated threshold function be h™ od ^ +1 ’ 2 \q). 
The transition map associated with each hybrid event is simply 
the identity map, q, ft, due to the continuity assumption. 
Our linearization of these hybrid dynamics towards an LTP 
approximation assumes that these transition times, t, zero 
crossings of h^q) and hi(q), maintain their periodicity and 
offsets within the period in close proximity of the limit-cycle, 
resulting in the following form of the nonlinear dynamics 


0=1 ( 6 ) 

. _ f fo(q,<t>,u) , if mod(f, T) G [0,t) 

1 \ /i(g,0,«) , if nrod(f, T) G [t,T) 

Assuming that the system given above has a limit cycle q(t) 
with a period T, linearization around q(t) yields the piecewise 
smooth LTP system 


x{t) 


Ao{t)x(t) + Bo(t)u(t), if mod {t,T) G [0,f) 
Ai(t)x(t) + Bi(t)u(t), if mod (t,T) G [t, T) 


where 
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u{t ) = 0 


q(t) = q(t) 
u(t ) — 0 


q(t) = q(t) 
u{t) — 0 


It is natural to assume that direct measurement of all x(t) 
may not be available or we may only measure a subset of x(t). 
Consequently, we also define a time-periodic output equation 
as in the form (9). 

Since system matrices Aft), Bft), Ci(t ) and Dft ) with 
i G {0,1} are time parametrized functions, the system has 
infinite parametric degrees of freedom, making parametric sys¬ 
tem identification challenging even when Harmonic Transfer 
Functions are used. At this point, we hypothesize that for 
hybrid systems, the variability within a chart is small compared 
to the change between charts and we approximate the LTP 
dynamics using a piecewise LTI approximation that preserves 


the LTP structure of the system. The LTP equations of motion 
then take the form 

i( t ) ~ / A ° X ^ + B ° U if mod (^ T ) G M) / 8) 

| Aix(t) + Biu(t), if mod(f,T) G [t,T) 

~ f C ox(t) + D 0 u(t ), if mod (t,T) G [0,f) 
y ' [ ) ~ 1 C\x(t) + D lU {t), if mod(f, T) G [t, T) 


The formulation above constitutes the basis of our framework 
for analyzing and identifying clock-driven legged locomotion 
models. 


III. Harmonic Transfer Functions 

A. Preliminaries and Background 

System identification studies on (stable) LTI systems rely 
on the fact that if the input is sinusoidal, then, at steady- 
state, the output will also be a sinusoidal signal (at the 
same frequency but with a possibly different magnitude and 
phase). This one-to-one mapping between input and output 
signals allows us to characterize the dynamics in terms of a 
frequency response function (FRF) also know as a Bode plot. 
Unfortunately, this approach does not readily transfer to LTP 
systems, which produce output spectra that include multiple 
(possibly infinite) harmonics of the input stimuli, each with 
possibly different magnitude and phase at steady state. 

One ad hoc way to mitigate this is to enforce a one-to-one 
mapping by neglecting higher harmonics [23]. However, this 
assumption may result in substantial inaccuracies particularly 
when the influence of higher harmonics on the response is 
expected to be significant. Motivated by this problem, Wereley 
[24] proposed a linear one-to-one mapping between the coef¬ 
ficients of an exponentially modulated periodic (EMP) signal 
at the input of LTP systems to the coefficients of an EMP 
signal at their output. This linear operator that maps the input 
harmonics to the output harmonics of an LTP system is called 
a Harmonic Transfer Function (HTF) [13]. 

In the following section, we review the derivation of 
harmonic transfer functions as presented in [13] and [25], using 
the principle of harmonic balance starting from the state space 
representation of (2). 

B. Theoretical Derivation of Harmonic Transfer Functions 

Recall that the system matrices in (2) are all T-periodic. 
Consequently, they can be represented by an infinite Fourier 
series with pumping frequency u p = 2i t/T. For the system 
matrix A(t), we have 

OO 

A(t)= £ A n e ju » nt . (10) 

n =—oo 

The matrices B(t), C(t) and D(t ) can be similarly decom¬ 
posed. In addition, we can also expand the state and output 
vectors since they are both EMP signals. Substituting these 
expansions into (2) and applying the principle of harmonic 
balance as explained in [13], we obtain the harmonic state 
space representation as 

sX = (A - N)X + BU 

y = cx + vu , 1! 1 ' 



where the doubly infinite vectors representing the harmonics 
of the state, control, and output signals are 


C. Estimation of Harmonic Transfer Functions via Frequency 
Domain System Identification 
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and the doubly infinite input modulation matrix is 

J\f := blockdiag{ jnco p I}, V?i € Z, 


( 12 ) 


(13) 


which modulates the input frequency to different harmonic 
frequencies. Details on the derivations can be found in [13], 

The T-periodic dynamics matrix, A(f ,), is expressed in 
terms of its complex Fourier coefficients, {A n \n £ Z}, as 
a doubly infinite block Toeplitz matrix. 
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with a similar definition for B(t) in terms of its Fourier 
coefficients represented by {B n \n £ Z}, C(t) in terms of 
{C„\n £ Z}, and D(t) in terms of {D n \n £ Z}. 

This collection of doubly infinite matrices is called the 
harmonic state space model (7~LSS) of the system given in (2). 
However, it will also be useful to determine an explicit input- 
output functional relationship between the Fourier coefficients 
of the harmonics of the input, {u n \n £ Z}, and those of 
the output, {y n \n £ Z}. This relationship is represented by 
the harmonic transfer function, G(s), which is also an infinite 
dimensional matrix of Fourier coefficients, satisfying 

y = GU. (15) 


Based on ('ll), G can be computed as 

G = C[sI-{A-N)]- 1 B + V , (16) 

as long as the inverse within this equation exists. 

There are, however, two problems associated with the 
harmonic transfer function as stated above. First, it is not 
clear whether the inverse of the doubly infinite matrix in 
the definition of the harmonic transfer function will always 
exist. This problem will be dealt with by an application of 
the Floquet Theorem. Second, the harmonic transfer function 
is a doubly infinite matrix operator, which cannot practically 
be implemented on a computer. This second problem will be 
mitigated by truncating the HTF in order to implement analysis 
on a computer. 

Note that the theoretical definition of harmonic transfer 
functions in [13], reviewed in this section, requires the state 
space representation of the system to be available. Our goal is 
to estimate this theoretically computed transfer function G(s) 
by using input-output data in the frequency domain without 
necessitating knowledge of internal system dynamics. 


In this section, we briefly explain the data-driven system 
identification method presented in [14]. 

In an LTP system, a sinusoidal input at a specific frequency 
generates a superposition of sinusoids at multiple (possibly 
an infinite number of) harmonics. Consequently, the system 
identification framework starts with tmncating the number of 
harmonic transfer functions G to be estimated. In the following 
examples, we consider only three frequencies in the output to 
clearly illustrate the approach of [14], 

Suppose that an LTP system consists of the superposition 
of three different harmonic transfer functions. Go, G_i and 
G i, each corresponding to a different frequency component of 
the output. The output can then be expressed as 

Y(jijj) = G 0 (ju])U (jut) + G-i(juj)U(juj + ju p ) 

+G 1 {jw)U(jw - ju p ). (17) 

In this new formulation, the n th transfer function is defined 
as the linear operator that maps the output at frequency io to an 
input at the same frequency, modulated with e? nUpt . However, 
a single input-output pair in each frequency will naturally 
not be sufficient to estimate harmonic transfer functions as 
in the case of LTI systems, since the identification problem 
will then be underdetermined. Therefore, either at least three 
independent inputs or additional constraints must be provided 
to enable a successful identification of these harmonic transfer 
functions. 

There are two key issues that need to be addressed before 
designing input signals for the identification process. First, 
we will require the use of at least as many variations on the 
input signal as the number of harmonic transfer functions to 
be estimated. This is accomplished in [14], which uses a single 
input sequence signal for system identification, constructed by 
concatenating phase shifted copies of a single waveform on the 
input evenly separated by delays within the system period. A 
complete characterization of system dynamics is possible with 
this method since different modes of the system were activated 
through the use of phase-shifted copies of a single waveform. 

The second issue is the need to excite all frequency compo¬ 
nents within the system by providing input signals with a suf¬ 
ficiently wide frequency spectrum. This can be accomplished 
through the use of chirp signals, whose frequency varies with 
time. The use of chirp input signals, combined with the idea of 
supplying multiple, phase-shifted input sequences allows us to 
obtain sufficiently rich input-output data to support the system 
identification process. 

Using this data with input-output pairs, one can estimate 
the harmonic transfer functions of the system, so that the error 
between actual and estimated outputs is minimized. Therefore, 
we can convert the identification problem to an optimization 
formulated as 

minimize J = (Y — U T G) 2 . ng) 

G 

However, note that Siddiqi [14] combines all phase-shifted 
signals in a single input. Hence, the problem is still underde¬ 
termined in the frequency domain, since a single input-output 









pair for a specific frequency will not be sufficient to identify 
three harmonic transfer functions. In order to address this 
problem, they consider additional constrains on the estimated 
harmonic transfer functions. First, they assume that transfer 
functions are smooth, which is reasonable for physical systems. 
This is enforced through a difference operator, D 2 , designed 
to compute the second derivative of a vector when multiplied 
from the left side. Details on the derivations for D 2 can be 
found in [14], 

The smoothness condition on transfer functions requires 
penalizing the curvature of individual transfer functions. There¬ 
fore, [14] modifies (18) to include a cost associated with the 
curvature, yielding a revised minimization problem formulated 
as 

minimize J = (Y — U T G) 2 + (aD 2 G) 2 , nm 
6 v J 

where a is a manually tuned constant weight for penalizing 
curvature. The solution of (19) can easily be obtained by 
differentiating J with respect to G, taking the form 

G = (U T U + aD 4 )- 1 U T Y , (20) 


where the rows of the matrix G (w) correspond to individual 
different harmonic transfer functions as 


GM 


' GiM ' 
GoM 
G-iM 


( 21 ) 


Fig. 2 illustrates the vertical leg model we focus on in this 
section. It consists of a mass attached to a leg with a spring- 
damper mechanism as well as a force transducer. Unlike the 
SLIP model, we assume that the toe is permanently affixed 
on the ground. Nevertheless, we recover the hybrid nature of 
locomotory gaits by assuming that the damper is turned on 
during a “stance phase” (lossy) and off during a “flight phase” 
(lossless). This construction recovers the hybrid nature of the 
dynamics, while allowing active input throughout the entire tra¬ 
jectory to support the generation of system identification data, 
as well as admitting theoretical computation of its harmonic 
transfer functions for a comparative investigation. 

We use the force transducer / in this system for two 
purposes. Firstly, active energy input to the system must be 
provided to maintain the limit cycle and compensate for energy 
losses due to the presence of damping. Second, it will be 
used as an exogenous input to the system to support the 
system identification process. Many physical legged platforms 
include similar active components in their legs to regulate 
their mechanical energy [26,27], Notwithstanding differences 
in how these actuators are incorporated into the system, they 
can all be used as the necessary exogenous inputs to perform 
system identification. A similar model was also investigated 
in [16] but using an additional nonlinear spring for energy 
regulation. 

The equations of motion for this simplified legged loco¬ 
motion model are given by 


Note that U T U and U T Y correspond to power spectral 
and cross spectral density functions, respectively. Therefore, 
(20) is analogous to estimating transfer functions in LTI 
systems, with an additional cost on curvature. 

IV. Simplified Legged Locomotion Model with 
Hybrid System Dynamics 

In this section, we describe a simple, vertically constrained 
spring-mass-damper system that possesses hybrid structural 
properties similar to the extensively studied Spring-Loaded 
Inverted Pendulum (SLIP) model for running behaviors. This 
will provide a simple example to illustrate the application of 
our system identification method to such systems. 


.. _ [~mg - cx - k(x - x 0 ) + f(t), if x > 0 
mX \—g — k(x — Xo) + f(t), otherwise. 

The lossy and lossless dynamics in (22) correspond to different 
charts in Fig. 1 and zero crossings of x represent threshold 
functions for both phases. 

Our illustrative examples use the parameters g = 9.81, k = 
200, c = 2, to = 1 and Xq = 0.2, chosen to be similar to the 
parameters of a vertical hopper platform in our laboratory [28]. 
As noted above, we choose the linear actuator input /(f) = 
/o(t) +u(t), consisting of a forcing term /o(f) to compensate 
for energy losses, and a chirp signal u[t) to introduce small 
periodic perturbations for system identification. 


A. System Dynamics 



Fig. 2. Simplified leg model. 


B. Theoretical Computation of Harmonic Transfer Functions 

The goal of this section is to compute harmonic transfer 
functions for our model around its limit cycle as outlined in 
Section III-B. 

We first assume that the forcing input / 0 (f) is appropriately 
chosen to induce an asymptotically stable limit cycle for this 
system. For example, our simple leg model achieves a stable 
limit cycle with /o(f) = cos(27rf). At this point, changing 
into error coordinates away from the limit cycle with / = 
x(t)—x{t), and substituting into (22), the equations of motion 
take the form 

-f-ct-kf, if £ + x(t) > 0 (23) 

kf, otherwise 

Due to the simplicity of the dynamics, this corresponds to 
a piecewise LTI system without necessitating any additional 
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approximations, taking the form 

J u(t), (24) 

where the hybrid nature of the system is captured by the flag 
s(£, f), with s = 1, when £ + x(t) > 0 and s = 0 otherwise. 


a 


0 

1 

a 
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—k 


a 


We now need to represent this piecewise LTI system as a 
linear time periodic system. However, even though the binary 
valued function s(£,t) can be considered time-periodic on the 
limit cycle itself, this is not the case for trajectories away 
from the limit cycle. To proceed, we hence assume that input 
induced perturbations are small, and that the binary valued 
function s(£,f) maintains its period and becomes strictly 
time dependent rather than state dependent, taking the form 
s(£, t) ~ s(t). We now can perform a Fourier series expansion 
on s(t) by treating it as a square wave with an offset to obtain 
a linear time periodic system in the form 


a' 


' 0 1 

'a' 

1 

'o' 
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—k —cs{t ) 

a 


1 


i(t), (25) 


y = t 1 o] 


Plugging these equations into the HTF framework described in 
Section III-B, yields analytic solutions to the harmonic transfer 
functions. We omit the details of this derivation due to space 
considerations, but use the resulting analytic solutions for the 
harmonic transfer functions up to rih = 10 to evaluate the 
output of our system identification method. 


C. Data-Driven Identification of Harmonic Transfer Functions 

In this section, we obtain harmonic transfer functions 
corresponding to the linearized dynamics of (25) by using 
input-output data without assuming prior knowledge of the 
state space model. Using / 0 (f) = cos(27rf) and u(t) = 0 for 
30 cycles without a perturbation, our example system stabilized 
to a limit cycle x(t) with a period T = Is. We use the 30 tft- 
period as the numerical limit cycle of the nonlinear system and 
subtract it from the trajectories of subsequent experiments to 
obtain the error function £i. 

In order to obtain input-output data for system identifica¬ 
tion, we apply an input signal consisting of nine subsequent 
30s long chirp signals, each with a linearly increasing fre¬ 
quency in the range (0, 7] Hz over its duration but with a 
different starting phase evenly distributed across the system’s 
period, T = Is. Each chirp signal has an amplitude of 0.004, 
chosen to be large enough to perturb system dynamics but 
small enough to keep the system close to the periodic orbit. 
A sample chirp signal with zero phase can be generated by 

u{t) = 0.004 sin(147rf 2 /30). (26) 


The resulting output is then subtracted from the numeri¬ 
cally measured limit cycle to obtain error trajectories for 
vertical position. The input signal and are then used as 
in Section III-C to estimate harmonic transfer functions for 
our system. Since our theoretical computations showed that 
responses beyond the third harmonic were very small, we only 
consider the fundamental harmonic and three harmonics on 
both sides for our experiments. 


-40 


-60 




Fig. 3. Estimation results for the fundamental harmonic. 


Fig. 3 illustrates the estimation performance of our al¬ 
gorithm for the magnitude and phase of the fundamental 
harmonic. Both graphs show that the application of the identifi¬ 
cation algorithm in [14] works well even for nonlinear periodic 
systems with hybrid dynamics. 

We also show our identification results for three harmonics 
in both the negative and positive sides in Fig. 4. Even though 
magnitudes for the harmonic transfer functions are small 
compared to the fundamental, the identification algorithm can 
provide accurate estimates for these transfer functions except 
in some narrow regions of G_ 2 and G 2 . The identification 
algorithm could not correctly estimate these two harmonics 
around 12—15 (rad/s). One possible reason for this discrepancy 
is the presence of strong responses in all harmonics around 
the same frequency except G _2 and G 2 , resulting in the 
inability of the identification algorithm to distinguish between 
the contributions from each harmonic absent knowledge of the 
internal system dynamics. Alternatively, these discrepancies 
may also be a result of the fact that hybrid transitions are not 
strictly time periodic (rather, they are state-dependent) which 
likely has effects on different frequencies and harmonics. We 
plan on investigating these issues further in the future. 

For a comparative analysis, we also present results from 
a parametric identification in order to show that further cor¬ 
rections on estimation results from a non-parametric method 
are possible. To this end, we fit the system parameters k and 
c in (25) by comparing root mean square error between theo¬ 
retically computed and estimated harmonic transfer functions 
Go, G_ 1 and Gi. We truncate the system response after the 
first harmonic in order to discard erroneous regions in higher 
harmonics. The resulting estimates were k = 200 for the 
spring constant and c = 2.12 for the damping coefficient, 
which closely coincide with the parameters used to generate 
the input-output data. As such, harmonic transfer functions 
obtained from parametric identification were found to closely 
match those obtained from theoretical computations as seen in 
Fig. 4. 

Motivated by these identification results, we plan to extend 
our work to the identification of the Spring-Loaded Inverted 
Pendulum (SLIP) model [3] and its extensions, widely used as 
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Fig. 4. Estimation results for the higher order harmonics. 


models of locomotory behaviors in the literature. Our future 
goal is to apply our system identification methods to our phys¬ 
ical monopod robot platform and to compare the identification 
performances with our previously verified analytical model 
[ 11 ]. 

V. Conclusion 

In this paper, we presented a system identification strat¬ 
egy to estimate input-output transfer functions for a simple 
hybrid spring mass damper system as a step towards data- 
driven models for legged locomotion. We first showed that 
a class of hybrid locomotion models can be approximated 
with a piecewise constant LTP systems in close proximity 
to their asymptotically stable limit-cycle. Our analysis and 
identification framework is based on the concept of harmonic 
transfer functions [24]. 

We first observed that the hybrid system dynamics associ¬ 
ated with this model exhibits piecewise LTI behavior around 
its periodic orbit. We then represented this behavior as a 
purely time periodic system around the limit cycle in order 
to utilize system identification techniques applicable to Linear 
Time Periodic systems. 

In order to provide a basis for comparison, we computed 
analytic expressions for harmonic transfer functions associated 
with the LTP approximation to our simplified hybrid model. In 
our theoretical analysis, we considered the system’s response 
up to the 10 th harmonic. We observed that there were no 
meaningful responses on both positive and negative sides after 
the third harmonic. Therefore, we decided to truncate the sys¬ 
tem response after the third harmonic during our identification 
studies. 

We then performed systematic simulation studies and iden¬ 
tified the harmonic transfer functions of the same model 
without knowledge of its internal dynamics. We used an 
input signal consisting of successive chirp signals, with phases 
evenly distributed across the system’s period, to obtain a full 
characterization of system dynamics for our frequency range 
of interest. Our studies showed that LTP system identification 
techniques can successfully be used to identify the transfer 


functions of nonlinear periodic models with hybrid system 
dynamics. 
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